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Dynamic Iterative Pursuit 

Dave Zachariah, Saikat Chatterjee and Magnus Jansson 

Abstract — For compressive sensing of dynamic sparse signals, we 
develop an iterative pursuit algorithm. A dynamic sparse signal process is 
characterized by varying sparsity patterns over time/space. For such sig- 
nals, the developed algorithm is able to incorporate sequential predictions, 
thereby providing better compressive sensing recovery performance, but 
not at the cost of high complexity. Through experimental evaluations, 
we observe that the new algorithm exhibits a graceful degradation 
at deteriorating signal conditions while capable of yielding substantial 
performance gains as conditions improve. 

I. Introduction 

Compressive Sensing (CS) |1| problems assume a sparse-signal 
model, undersampled by a linear measurement process. The algo- 
rithms for CS can be separated into three broad classes: convex 
relaxation, Bayesian inference, and iterative pursuit (IP). For large- 
dimensional CS signal-reconstruction, IP algorithms offer compu- 
tationally efficient solutions. Examples of such IP algorithms are 
orthogonal matching pursuit (OMP) [2|, subspace pursuit (SP) t3l 
and several variants of them [4|[5 1[6 |[7 ']. The methodology of such IP 
algorithms is to detect and reconstruct the non-zero, or 'active', signal 
coefficients in a least-squares framework. These algorithms may use 
some prior information, such as the maximum allowable cardinality 
of the 'support set'. The support set is defined as the set of active 
signal coordinates of the underlying sparse signal. In general, the IP 
algorithms work with a single snapshot of the measurements. In this 
paper, we are interested in generalizing the iterative pursuit approach 
so as to use more prior information. Such prior information may for 
instance be available in dynamically evolving sparse processes with 
temporal/spectral/spatial correlations, as in the sparse signal scenarios 
of magnetic resonance imaging (MRI) (8), (9J, spectrum sensing 1 10 1 
and direction of arrival estimation Bill . 

Incorporation of prior information is a recent trend in CS. In 1121 . 
the overall methodology is sequential and can be seen as a two- 
step approach: (1) support-set detection of the sparse signal, and 
(2) reduced-order recovery using prior information on the detected 
support set. For a reasonable detection of support set, 0121 uses 
convex relaxation algorithms. Then, a standard Kalman filter (KF) 
is employed to use prior information for sequential signal recovery. 
Without explicit support set detection, 1131 uses KF to estimate 
the entire signal and enforces sparsity by imposing an approximate 
norm constraint. However, the work of |13] validates their algo- 
rithm for a signal with a static sparsity pattern (i.e. an unknown 
pattern that does not evolve over time). Similarly, 11141 considers 
scenarios with static sparsity patterns and solves the reconstruction 
of a temporally evolving sparse signal with multiple measurement 
vectors in a batch Bayesian learning framework with unknown model 
parameters. Iterative pursuit algorithms that can use prior information 
to recover dynamic sparse signals are, however, largely unexplored. 
One exception is 1151 which uses a maximum aposteriori criterion 
to modify SP for Gaussian processes. Their signal model, however, 
does not allow explicit modeling of the temporal correlation of the 
sparsity pattern. 

In this paper, we consider a signal model with dynamically 
evolving sparsity pattern. In other words, we consider that the signal 
sparsity pattern varies over time/space at any rate (i.e. from a 
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slowly varying case to a rapidly varying case). We then develop a 
predictive orthogonal matching pursuit algorithm that can incorporate 
prior information in a stochastic framework, using the signal to 
prediction error ratio as a statistic. Thereby recovery performance 
can be improved while maintaining the complexity advantage of IP 
algorithms. This generalizes the linear minimum mean square error 
(MMSE) approach taken in Gaussian-based matching pursuit 1161 . 

We also develop a robust detection strategy for finding the support 
set elements of a dynamic sparse signal. Compared to standard 
correlation-based successive detection in existing iterative pursuit 
algorithms, such as OMP, this detection strategy is found to be more 
robust to erratic changes in the sparsity pattern. 

Finally, the algorithm is integrated into a recursive Kalman-filter 
framework in which the sparse process is predicted as a superpo- 
sition of state transitions. The new IP algorithm using sequential 
predictions is referred to as dynamic iterative pursuit (DIP). Through 
experimental simulations we show that the new algorithm provides a 
graceful degradation at higher measurement noise levels and/or lower 
measurement signal dimensions, while capable of yielding substantial 
gains at more favorable signal conditions. 

Notation: ||x||o denotes lo 'norm', i.e. the number of non-zero 
coefficients of the vector x. A©B is the direct sum of matrices. |S| 
and S c are the cardinality and complement of set S, respectively. 
denotes the empty set. (•)* is the Hermitian transpose operator. A' 
the Moore-Penrose pseudoinverse of matrix A. C 1 ' 2 denotes a matrix 
square root of a positive definite matrix C, and C*'' 2 is its Hermitian 
transpose. A^ x ,j\ denotes a submatrix of A with elements from row 
and column indices listed in ordered sets I and J . Similarly, the 
column vector xpq contains the elements of x with indices from set 
X. 

II. Signal model 
We consider a standard CS measurement setup, 

y t = Hx t +n t G<C M , (1) 

where xt G is the sparse state vector to be estimated and the 
tit is zero-mean Gaussian, ElritnjL;] = R t <5(7). The sensing matrix 
H = [hi ■ • ■ hjv] G C MxJV , where M < N. Both H and R t 
are given. Without loss of generality, we assume ||h.;||2 = 1. 

A. Process model 

Let the 'support set' I X]t C {1, ...,N} represent the sparsity 
pattern of x t G C . It will be assumed that \I x ,t\ = ||xt||o < K max , 
where K mllx < M. Let Xji denote the state transition probability 
j — > i of the 'active' signal coordinate j. Then the probabilities 
determine the transition I Xt t — > I x ,t+i, as will be illustrated below. 

The transition of an active signal coordinate j — > i is modeled as 
an autoregressive (AR) process, 

Xi, t+ i = ctijXj, t + Wi,t, (2) 

where Xj,t denotes the jth component of x t , Wi, t is the associated 
innovation and the AR coefficient \ctij\ < 1. This model extends 
the scenario considered in [ 14| where the transition probabilities are 
degenerate Xji = S(i — j), resulting in a static sparsity pattern I Xt t = 
i»,Vt. 

The sparse-signal process can be written compactly as a linear 
state-space model with random transition matrices At and B t , 

x f+ i = A t xt + B t w t , (3) 

where wt is zero-mean Gaussian, E[w t W(_ ; ] = QS(l) G C NxN 
and Q = diag(cr 2 , . . . , o%). The non-zero elements of A t G C NxN 
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Fig. 1. Example of evolving sparsity pattern with (N,K,T) 
(200, 10, 200) and transition probabilities (|4j. 



are aij,t = otij for all j G I x ,t and i G I x ,t+i- Similarly, the non- 
zero elements of the diagonal matrix B f G 
all i (E I x ,t+i. The model parameters ctij, X 
be known. 



l N * N are 6», t = 1 for 
j i and Q are assumed to 



B. Examples 

Use of the transition probabilities Xji along with signal model l(2j 
enables compact modeling of dynamically evolving sparsity patterns. 
The potential applications include MRI, spectrum sensing, direction 
of arrival estimation, frequency tracking etc. 

As an initial example, consider a slowly varying sparsity pattern 
I x ,t over T snapshots, following 

'0.90 i = j 

0.05 i=j±l, ifjg{l,N} 

Xji = { 0.10 i = j + l, if j = 1 . (4) 

0.10 i=j-l, if j = N 

\i-j\>l 

A realization of this process is illustrated in Figure [TJ This choice is 
intended to model the strong temporal correlation of sparse signals 
exhibited in e.g. MRI. 

Next, consider a simpler parameterization, 



X, 



1 



1 = 3 



(5) 



where v G [0, 1] is a mixture factor. This is intended to model 
more erratically evolving patterns in e.g. frequency-hopping radio 
frequency (RF) signals. Examples of resulting sparsity patterns are 
shown in Figure [2] where we consider v = 0.01 and 0.5. It can be 
seen that the evolution of the sparsity pattern becomes more erratic as 
v increases. In the above examples we have ensured that the sparsity 
level is constant, K = 10. 

III. Dynamic Iterative Pursuit 

We approach the dynamic estimation problem by first developing 
an iterative pursuit algorithm that can incorporate prior information 
in the form of a prediction of xt. A support set detection strategy is 
proposed using the signal to prediction error ratio. Next we develop 
a recursive algorithm based on the Kalman-filter framework. We 
propose predicting the sparse process as the superposition of all state 
transitions of the signal coefficients. 

'The two cases j = 1 or j = N are necessary for the edge states. 
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Fig. 2. Examples of evolving sparsity patterns with (TV, K, T) = 
(200, 10, 100) and transition probabilities j5] with (a) v = 0.01 and (b) 
v = 0.5. 



A. Incorporation of prior information 

Given the constraint on the support set, |/a,| < JTmax, the brute 
force least-squares solution would be to enumerate all combinations 
of possible support sets. For each set, I C {1, . . . ,N}, the signal 
coefficients are reconstructed by a least-squares criterion and a mea- 
surement residual is computed, r = y-Hi. nxm, The reconstruction 
with minimum residual norm is then chosen as the solution. However, 
with at least one active coefficient there are YIkSl (jD possible 
support sets I to enumerate, which is clearly intractable. 

Several iterative pursuit algorithms solve the estimation problem 
by a sequential detection of the support set and reconstruction of the 
corresponding signal coefficients. We will use OMP to illustrate the 
essential components of this sequential strategy. 

OMP takes a support set / as its starting point. Reconstructed 
signal coefficients, Xj, j G /, are cancelled from the observation y 
to form the residual r = y — Hr. nxrji. Initially 1 = 0. Under the 
hypothesis of a remaining active coefficient Xi, i J, the residual 
signal model is 

r = hiXi + + n, 

where £j = Xj — Xj are estimation errors. OMP detects the active 
coefficient by using a matched filter. The matched filter employs the 
strategy of estimating Xi using a least-squares criterion, Xi = hjr = 



(6) 



h*r. The index i £ I corresponding to maximum energy is 
added to /. Finally the coefficients corresponding to / are estimated 
jointly based on a least-squares criterion, solving 



xrn = argmin 



H 



[.,ijx[i]|| 2 = H [.,/]y 



The residual r is updated and the process is repeated until the residual 
norm no longer decreases or when \I\ reaches the limit K m!a . For 
sake of clarity OMP is summarized in Algorithm [TJ where k denotes 
the iteration index. 

Using a stochastic framework, we now extend the estimation 
strategy to a scenario in which a prediction x~ = x + e is given, 
where e ~ A/"(0,P~) and error covariance matrix P~ is known. 
Then the signal to prediction error ratio, 

^S^EUlCoo), (7) 



E[|e,| 2 ] 

quantifies the certainty that i belongs to the support set. We propose 
to use pi for selecting indices to be added to I. The ratio pi is 
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Algorithm 1 : Orthogonal Matching Pursuit (OMP) 



Algorithm 2 : Predictive Orthogonal Matching Pursuit (PrOMP) 
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Given: y and H 

Set k — 0, r = y and I = 

repeat 

k := Jfc + 1 

ik = argmax i6J c |h*r fc _i| 
I :=IUi k 



Hi] 



H 



7 ,y; *[/c] = o 



r k = y - H[.,j]X[z] 
until (||r fc || 2 > ||r fc _i|| 2 ) or (k > K m!tK ) 
Output: x and / 



successively updated by conditioning the expectations on the residual, 
under the hypothesis with signal model ([6). Then E[\xi | 2 ] = \pi\ r | 2 + 
<7 2 | r where the conditional mean //j| r is given by the MMSE-estimator 
and (jf| r by its error variance. The prior of Xi is the prediction x~ . 

For tractability the estimation errors £ j are assumed to be Gaussian 
and their correlations negligible so that the MMSE-estimator gives 



H'i | r 
2 



^ + Si (r - h t x- 
■ (1 - gihOpf, 



(8) 



where p 4 is the ith diagonal element of P 1171 , The gain (row) 
vector g* and covariance matrix D are, 



g* = (— +h*D- 1 h l N ) h*D 1 



Pi 

D = ^Vfh,h* +R, 

ie/ 



(9) 



where <jj is the variance of £j and R is the covariance matrix of 
n. As the support set I successively grows, the inverse D 1 can be 
updated efficiently using the Sherman-Morrison formula, as shown 
below. 

To sum up, the signal to prediction error ratio is given by 



Pi 



I i2 i 2 
lMi|r| + °~i\r 



(10) 



and approximated using ([8}. For the maximum pi, i £ I is added 
to I. Finally, signal coefficients are jointly re-estimated, solving a 
weighted least-squares problem 



X[j] = argmin 



X[7] 



(11) 



where S = PJ^ „. This is the linear MMSE estimator provided I is 
the correct support set. The residual r is updated and the process is 
repeated as above. The resulting algorithm is referred to as 'Predictive 
OMP' (PrOMP) and is summarized in Algorithm [2] 

The function mmse-rec solves Jilt and can be computed by a 
measurement update of form: 



XM = X[I] +K(y-H hI]X[/] J , 

where K = (S _1 + H^R^H^^Hj^R" 1 . The function 
update-cov updates the inverse covariance matrix for the added 
reconstructed coefficient Xi, 



D 1 :=D 1 I 



h.KD- 1 
ar 2 + hi-D-ihi 



where of is the corresponding diagonal element of the posterior error 
covariance matrix (S _1 +Hf J ]R~ 1 H[.,/]) _1 1171 . 



Given: y, H, R , x and P 

Set k = 0, r = y, / = and D" 1 = R -1 

repeat 

k ~ k + 1 

Compute pi using d 10b and ^ 

ik = argmax i6/ c pt 
I:=IUi k 

X[j] =mmse-rec(y,H,R _1 ,x _ ,P _ ,7): 
r k = y - H[.,j]X[/] 

D" 1 =update-cov(D- 1 ,H,R- 1 ,p-, 
until (||r fe || 2 > ||r fc _i|| 2 ) or (fc > K mm ) 
Output: x and I 



*[/<=] 
•f,ifc) 



B. Robust support-set based strategy 

The strategy described above performs a successive cancellation 
of reconstructed signal coefficients. The performance is therefore 
crucially dependent on detecting an active coefficient individually 
at each stage, which is a 'hard' decision. But the hypothesis of one 
remaining active coefficient at each stage induces a risk of irreversible 
detection errors. This increases with more erratically evolving sparsity 
patterns, since the process is harder to predict. 

The signal to prediction error ratio pi, however, provides a statistic 
that can be viewed as 'soft information'. In order to increase 
robustness to detection errors we propose to use pi for selecting 
the t most likely remaining coefficients. Let us denote the set of I 
most likely indices by L. It is joined with the existing set / to form 
a hypothesized support set / = / U L. This set is used to reconstruct 
Xrfi and the coefficient Xi, i € L with maximum magnitude is added 
to the support set / at each stage. Here I = max(0, K mdx ~ |/|), which 
prevents overfitting beyond the prior knowledge of the sparsity level. 

Algorithm [3] describes this alternative detection strategy, based on 
a hypothesized support set. The concerned scheme is called 'robust 
predictive OMP' (rPrOMP). 

Algorithm 3 : Robust predictive OMP (rPrOMP) 



R 



Given: y, H, R , x and P 
Set k = 0, r = y t , I = and D _1 
repeat 

k := k+1 

Compute pi using d 10b and l[8} 

£ = max(0, K mi x - \I\) 

L — {indices of £ largest pi G I c } 

I = IUL 

Xrj, =mmse-rec(y, H, R _1 , x~, P~ , 7) 
i k = argmax, €i \x l \ 
I := I Uik 

X[j] =mmse-rec(y,H,R~ 1 ,x~,P~,7); X[jej = 
v k = y - H^/jx^] 

D _1 =update-cov(D _1 , H, R 1 , P" , I, i k ) 
until (||r fe || 2 > ||r fc _i|| 2 ) or (k > K max ) 
Output: x and I 



C. Prediction of dynamic sparse signals 

Suppose a snapshot yt has been observed and a prediction x<T is 
given along with P^ • Let x t denote the estimated sparse state vector 
after the application of a predictive greedy pursuit algorithm (either 
PrOMP or rPrOMP), and / its support set. Then the updated error co- 
variance matrix Pt is computed block-wise corresponding to the set I 
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and its complement I c . First, P^jj.t = (S^ 1 + H*. 7 jR7 1 H[.,/]) 
is the posterior error covariance, where St = PJ} j, t 1171 . Sec- 
ond, the uncertainty of the inactive coefficients is preserved by 
P[/ c ,i c ],t = P[}c jc] t - Finally, in line with the MMSE reconstruction 
jilt , the cross-correlations are set as Prr /ci t = and Prfc n t = 0. 

We propose predicting Xt+i from Xt as a superposition of all 
possible transitions, 



or written compactly, x t 



N 
3 = 1 

= Fxt, where /y 



A, 



(12) 



The 



prediction error covariance matrix is then approximated by the 



equation, P 



t+i 



FP t F* + Q. 



Putting these blocks together we develop a Kalman-filter based 
algorithm for recovery of sparse processes in Algorithm which we 
call dynamic iterative pursuit (DIP). In DIP we use predictive OMP 
(PrOMP). If robust predictive OMP (rPrOMP) is used instead, the 
algorithm can be referred to as 'rDIP'. 

Algorithm 4 : Dynamic Iterative Pursuit (DIP) 



for t = 0, . . . do 

%Measurement update 

[xt,J] -PrOMP(yt,H,R t -\x t -,P t -) 



P[i<=,i<=],t = 
%Prediction 

K+i = F *t 

FP t F* + Q 



= 0; P [7 



= 



end for 



In the experiments, two signal parameters are varied; (a) the signal- 
to-measurement noise ratio, 



SMNR- E E*I« 



EEtl|n t ||| 
while fixing E[||xt||l] = 1 so that a\ — 



(13) 

MX SMNR ' and the 



fraction of measurements k = M/N. 

For a performance measure we use the signal-to-reconstruction 
error ratio, defined as 



SRER = 



E[EJ|x t ||i] 



EE, 



Mil 



(14) 



which is the inverse of the normalized MSE. Note that SRER = dB, 
i.e. no reconstruction gain, is equivalent to using x t = 0. 

B. Algorithm initialization 

For the predictive algorithms — DIP, rDIP and genie-aided KF — 
we use the mean and variance of an autoregressive process as initial 



values, x n 



and Pn 



<t x In where a x 



In these 



v Q — u aiiu x q u x ±n wucic u x — 1 _ a ' i . 

algorithms we set K m!a — K for consistent comparisons, although 
strict equality is not a requirement. 

Here we mention that BPDN 1 1 ] solves 



xt = argmin ||xt||i subject to ||yt 



Hx t || 2 < s, 



where the slack parameter e is determined by the measurement noise 
power, as 



y<r2(M + (2v / 2M)), 



following 1181 . 1161 . Note that BPDN does not provide a A"-element 
solution. It is also unable to use prediction. The code for BPDN is 
taken from the Zi-magic toolbox. 



IV. Experiments and Results 

In this section we evaluate DIP with respect to static OMP, SP 
and convex relaxation based basis pursuit denoising (BPDN) |T) 
algorithms. We also show the performance of a 'genie-aided' Kalman 
filter (KF) which provides a bound for MMSE-based reconstruction 
of linear processes. The genie-aided approach is given the sparsity 
pattern a priori, but does not know the active signal coefficients. 
Finally, the robustness properties of rDIP are compared with DIP 
for erratically evolving sparsity patterns. The results are shown using 
Monte Carlo simulations, averaged over 100 runs. 

A. Signal generation and performance measure 

Using a typical setup we consider a sparse process with the 
parameters N = 200, K = 10 and number of snapshots T = 100, 
with oscillating coefficients according to an AR-model as in {2]l with 
Oiij — a = —0.8, and Q = cr^Ijv. The sparsity pattern transitions, 
Ix,t — > Ix.t+i, are determined by transition probabilities Xji which 
are set in the experiments. 

The transition of each active state j G I x .± is generated by a first- 
order Markov chain with Xji. If two states in I Xi t happen to transition 
into one, a new state is randomly assigned to I x ,t+i, to ensure that 
the sparsity level is constant in the experiment. 

The entries of the sensing matrix H are set by random drawing 
from a Gaussian distribution A/"(0, 1) followed by unit-norm col- 
umn scaling. The measurement noise covariance matrix has form 
Rt = ofjljvf. Process and measurement noise are generated as 
w t ~ N(0, Q) and n t ~ Af{0, R t ), respectively. 



C. Results 

For all experiments we ran 100 Monte Carlo simulations, where a 
new realization of {xt, yt}t=i and H was generated for each ran. 

In the first experiment we consider a slowly varying sparsity pattern 
I Xi t following the transition probabilities Xji in @. Figure [3] shows 
how the algorithms perform with varying measurement noise power 
at a fixed fraction of measurements k = 0.25. DIP overtakes static 
BPDN at lower SMNR levels, while exhibiting a similar graceful 
degradation. The static OMP and SP do not take into account the 
measurement noise and hence continue to degrade. For instance, DIP 
reaches the cut-off point of dB reconstruction gain at an SMNR 
level that is approximately 5 dB lower than the static OMP. Figure [4] 
shows how the improvements persist for varying k at a fixed SMNR 
= 10 dB. In this scenario rDIP exhibits similar performance as DIP. 
Taking static OMP as the baseline algorithm, the improvement of 
predictive iterative pursuit is illustrated in Figure |5] The minimum 
SRER advantage is about 2 dB, and increases substantially with rising 
SMNR. 

Next, we consider an unknown but static sparsity pattern I Xtt = I x , 
generated by degenerate transition probabilities Xji — S(i — j), 
and compare DIP with a 'genie-aided' KF. The latter filters the 
coefficients of a known support set I X) t, and therefore provides 
an upper bound on the performance of sequential estimation. The 
bound is not necessarily tight since only \I x ,t\ < K is given in the 
problem. As SMNR increases, DIP rapidly approaches the bound 
while OMP saturates for k = 0.25, illustrated in Figure [6] At SMNR 
= 20 dB, OMP and DIP are about 10 and 2 dB from the upper limit, 
respectively. Again, rDIP performs similarly to DIP. 
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Fig. 3. Comparison of different methods where we show SRER versus 
SMNR at re = 0.25 and set transition probabilities according to (@). 
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Fig. 5. SRER improvement over OMP versus SMNR at re = 0.25. Transition 
probabilities according to @. 




Fig. 4. SRER versus re = M/N at SMNR = 10 dB. Transition probabilities 
according to l(4j. 



Finally, we consider an erratically evolving sparsity pattern I x ,t, 
with transition probabilities Xji set according to l[5}, in order to 
compare the robustness of rDIP with DIP. Figure [7] shows how perfor- 
mance is affected as the mixture factor v increases. DIP converges to 
OMP from above; rDIP provides near equivalent performance to DIP 
at first but shows a more graceful degradation. At the extreme, when 
all transitions are equiprobable, rDIP is still capable of yielding above 
+2.5 dB gain over OMP. This validates the robustness considerations 
behind its design. 

Reproducible results: MATLAB code for the algorithms can be 
downloaded from http://sites.google.com/site/saikatchatt/softwares 
The codes produce the results in Figures [3] and UJ 

V. Conclusions 

We have developed a new iterative pursuit algorithm that uses se- 
quential predictions for dynamic compressive sensing, which we call 
dynamic iterative pursuit. It incorporates prior statistical information 
using linear MMSE reconstruction and the signal to prediction error 
as a statistic. The algorithm was experimentally tested on a sparse 
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Fig. 6. SRER versus SMNR at re = 0.25. Degenerate transition probabilities 
(fixed sparsity pattern). 



signal with oscillating coefficients and evolving sparsity pattern. The 
results show that the algorithms exhibit graceful degradation at low 
SMNR regions while capable of yielding substantial performance 
gains as the SMNR level increases. 
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